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ABSTRACT 

This  paper  discusses  fundamental  issues  related  to  the  integration  of  computer  aided  design  and 
analysis  (I-CAD-A)  by  introducing  a  new  class  of  ideal  compliant  joints  that  account  for  the 
distributed  inertia  and  elasticity.  The  absolute  nodal  coordinate  formulation  (ANCF)  degrees  of 
freedom  are  used  in  order  to  capture  modes  of  deformation  that  cannot  be  captured  using 
existing  formulations.  The  ideal  compliant  joints  developed  can  be  formulated,  for  the  most 
part,  using  linear  algebraic  equations,  allowing  for  the  elimination  of  the  dependent  variables  at 
a  preprocessing  stage,  thereby  significantly  reducing  the  problem  dimension  and  array  storage 
needed.  Furthermore,  the  constraint  equations  are  automatically  satisfied  at  the  position, 
velocity,  and  acceleration  levels.  When  using  the  proposed  approach  to  model  large  scale  chain 
systems,  differences  in  computational  efficiency  between  the  augmented  formulation  and  the 
recursive  methods  are  eliminated,  and  the  CPU  times  resulting  from  the  use  of  the  two 
formulations  become  similar  regardless  of  the  complexity  of  the  system.  The  elimination  of  the 
joint  constraint  equations  and  the  associated  dependent  variables  also  contribute  to  the  solution 
of  a  fundamental  singularity  problem  encountered  in  the  analysis  of  closed  loop  chains  and 
mechanisms  by  eliminating  the  need  to  repeatedly  change  the  chain  or  mechanism  independent 
coordinates.  It  is  shown  that  the  concept  of  the  knot  multiplicity  used  in  computational 
geometry  methods,  such  as  B-spline  and  NURBS  (Non-Uniform  Rational  B-Spline),  to  control 
the  degree  of  continuity  at  the  breakpoints  is  not  suited  for  the  formulation  of  many  ideal 
compliant  joints.  As  explained  in  this  paper,  this  issue  is  closely  related  to  the  inability  of  B- 
spline  and  NURBS  to  model  structural  discontinuities.  Another  contribution  of  this  paper  is 
demonstrating  that  large  deformation  ANCF  finite  elements  can  be  effective,  in  some  MBS 
application,  in  solving  small  deformation  problems.  This  is  demonstrated  using  a  heavily 
constrained  tracked  vehicle  with  flexible  link  chains.  Without  using  the  proposed  approach, 
modeling  such  a  complex  system  with  flexible  links  can  be  very  challenging.  The  analysis 
presented  in  this  paper  also  demonstrates  that  adding  significant  model  details  does  not 
necessarily  imply  increasing  the  complexity  of  the  MBS  algorithm. 

Keywords:  Compliant  joints;  tracked  vehicle;  flexible  multibody  systems;  finite  element; 
ANCF;  closed  loop  chains;  I-CAD-A 
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1.  INTRODUCTION 

Complex  multibody  systems  (MBS)  consist  of  many  components  connected  by  mechanical 
joints.  These  joints  can  be  modeled  using  different  formulations  that  lead  to  different  dynamic 
models.  One  method  for  representing  mechanical  joints  in  computational  MBS  algorithms  is  the 
ideal  joint  formulation,  in  which  a  set  of  nonlinear  algebraic  constraint  equations  are  used  to 
define  the  body  connectivity  conditions.  This  approach,  when  used  with  existing  MBS  rigid  body 
formulations  or  flexible  body  formulations  based  on  the  floating  frame  of  reference  (FFR) 
formulation,  does  not  account  for  the  joint  flexibility.  Furthermore,  in  these  existing  MBS 
formulations,  the  joint  constraints  introduce  nonlinearities  both  in  the  inertia  and  constraint 
forces.  Another  problem  that  must  be  addressed  when  constrained  MBS  dynamics  algorithms  are 
developed  is  the  violation  of  the  joint  kinematic  constraints  that  must  be  satisfied  at  the  position, 
velocity,  and  acceleration  levels.  This  paper  proposes  a  new  procedure  to  address  the  following 
fundamental  issues: 

1.  The  paper  develops  a  new  class  of  ideal  joints  that  account  for  the  distributed  inertia  and 
elasticity  at  the  joint  definition  point.  Additional  degrees  of  freedom  are  introduced  in 
order  to  capture  modes  of  deformation  that  cannot  be  captured  using  existing 
formulations  that  employ  rigid  triads  or  unit  vectors  to  formulate  the  joint  equations. 

2.  The  ideal  compliant  joints  developed  in  this  paper,  while  more  general  and  capture  more 
modes  of  displacements,  can  be,  for  the  most  part,  formulated  using  linear  algebraic 
equations.  These  linear  connectivity  conditions  allow  for  the  elimination  of  the  dependent 
variables  at  a  preprocessing  stage,  thereby  significantly  reducing  the  problem  dimension 
and  array  storage  needed. 
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3.  Because  linear  connectivity  conditions  are  used  to  eliminate  the  dependent  variables  at  a 
preprocessing  stage,  the  constraint  equations  are  automatically  satisfied  at  the  position, 
velocity,  and  acceleration  levels.  Therefore,  when  using  the  proposed  approach  to  model 
large  scale  chain  systems,  differences  in  computational  efficiency  between  the  augmented 
formulation  and  the  recursive  methods  are  eliminated,  and  the  CPU  times  resulting  from 
the  use  of  the  two  formulations  become  similar  regardless  of  the  complexity  of  the 
system. 

4.  Furthermore,  the  elimination  of  the  joint  constraint  equations  and  the  associated 
dependent  variables  contribute  to  the  solution  of  a  fundamental  singularity  problem 
encountered  in  the  analysis  of  closed  loop  chains  and  mechanisms.  It  is  known  that  in  the 
case  of  closed  loop  chains  and  mechanisms,  the  set  of  independent  coordinates  must  be 
repeatedly  changed  in  order  to  avoid  singularity  of  the  constraint  Jacobian  matrix.  This 
problem  can  be  avoided  by  using  the  new  ANCF  finite  element  meshes  that  allow  for  the 
elimination  of  the  constraint  equations  and  the  dependent  variables  at  a  preprocessing 
stage. 

5.  Using  the  concept  of  the  ideal  compliant  joints,  the  fundamental  problems  that  can  be 
encountered  when  using  computational  geometry  methods,  such  as  B -spline  and  NURBS, 
as  analysis  tools  are  explained.  The  concept  of  the  knot  multiplicity  used  in  computational 
geometry  methods  does  not  lend  itself  easily  to  the  formulation  of  the  ideal  compliant 
joints  that  may  require  different  degrees  of  continuity  for  different  variables  or  elements 
of  the  gradient  vectors  or  may  require  constraints  on  the  gradient  vectors  along  directions 
that  differ  from  the  B -spline  and  NURBS  coordinate  lines.  This  issue  is  closely  related  to 
the  inability  of  B-spline  and  NURBS  to  model  structural  discontinuities. 
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6.  While  small  deformation  MBS  problems  can  be  effectively  solved  using  the  FFR 
approach  that  allows  for  filtering  out  high  frequency  contents,  there  are  small 
deformation  problems  that  can  be  solved  more  efficiently  using  large  deformation 
formulation  such  as  the  absolute  nodal  coordinate  formulation  (ANCF).  In  these  small 
deformation  problems,  the  inclusion  of  high  frequency  modes  may  be  necessary.  An 
example  of  these  problems  is  a  flexible-link  chain  of  a  tracked  vehicle  where  each  link  is 
modeled  using  one  finite  element  only.  The  use  of  ANCF  finite  elements  can  be 
advantageous  as  compared  to  the  use  of  the  FFR  formulation.  This  important  issue  is 
addressed  in  this  paper. 

7.  The  necessity  of  using  the  new  concepts  developed  in  this  paper  in  modeling  complex 
systems  is  demonstrated  using  a  heavily  constrained  tracked  vehicle  with  flexible-link 
chains.  Without  using  the  proposed  approach  modeling  such  a  complex  system  with 
flexible  links  can  be  very  challenging. 

The  analysis  presented  in  this  paper  demonstrates  that  adding  significant  model  details  does 
not  necessarily  imply  increasing  the  complexity  of  the  MBS  algorithm.  In  fact,  including 
significant  details  such  as  body  flexibility  and  joint  deformations  can  be  accomplished  using 
MBS  algorithms  that  are  much  simpler  than  the  algorithms  currently  used  for  rigid  body  system 
analysis.  This  can  be  achieved  using  the  ANCF  finite  elements  that  lead  to  new  meshes  that  have 
linear  connectivity  and  constant  inertia.  By  using  the  geometry  concepts  that  allow  for  successful 
integration  of  computer  aided  design  and  analysis  (I-CAD-A)  and  using  the  new  ANCF  meshes 
to  model  heavily  constrained  systems  such  as  long  chains,  the  distinction  between  the  augmented 
and  recursive  methods  can  be  eliminated. 
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2.  BACKGROUND 

Several  FE/MBS  ,  formulations  have  been  developed  for  the  deformation  analysis  of  mechanical 
systems.  Some  of  these  formulations  are  suited  for  the  analysis  of  large  rotation  small 
deformation  problems,  while  the  others  can  be  used  in  the  analysis  of  large  rotation  and  large 
deformation  problems.  An  example  of  the  large  rotation  and  small  deformation  formulation  is  the 
FE/FFR  formulation  which  employs  mixed  set  of  reference  and  elastic  coordinates  (Shabana, 
2014).  This  formulation  leads  to  a  highly  nonlinear  inertia  matrix  that  has  a  strong  coupling 
between  the  reference  and  deformation  coordinates.  In  addition  to  the  high  nonlinearity  of  the 
inertia  forces,  the  joint  constraints  equations  developed  using  the  FFR  formulation  are  also 
highly  nonlinear  because  of  the  use  of  orientation  parameters  (Korkealaakso  et  al,  2009). 
Furthermore,  rigid  triads  and  unit  vectors  are  often  used  in  the  formulation  of  the  joints  in  the 
rigid  body  and  FFR  algorithms,  and  as  a  result,  these  ideal  joint  formulations  do  not  capture 
modes  of  deformation. 

While  the  FFR  formulation  will  remain  an  effective  approach  for  modeling  small 
deformation  problems  because  it  allows  for  filtering  out  high  frequency  modes,  the  more  general 
large  rotation  large  deformation  finite  element  absolute  nodal  coordinate  formulation  (ANCF) 
can  be  effective  in  solving  a  class  of  small  deformation  problems  (Shabana  et  al,  2012;  Hamed  et 
al,  2011;  Shabana,  2010,  2012;  Dimtrochenko  and  Pogorelov,  2003;  Abbas  et  al,  2010;  Dufva  et 
al,  2005,  2007;  Schwab  and  Meijaard,  2010;  Shabana  and  Mikkola,  2003;  Tian  et  al,  2009; 
Yakoub  and  Shabana,  2001;  Omar  and  Shabana,  2001).  ANCF  finite  elements  can  be  used  to 
model  both  structural  and  nonstructural  discontinuities  in  a  straightforward  manner  (Shabana  et 
al,  2012;  Hamed  et  al,  2011).  Using  ANCF  finite  elements,  the  absolute  position  vectors  and 
their  gradients  are  considered  nodal  coordinates.  The  global  position  of  an  arbitrary  point  on  an 
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element  k  of  body  i  can  be  defined  in  a  global  coordinate  system  as  r‘k  =  S'k(x,k,  y'k ,  z,lk)e,k  (t) , 
where  S ,k  (x'k ,  y'k ,  z'k )  is  the  space-dependent  matrix  that  defines  the  element  shape  functions, 

and  e'A  is  time-dependent  vector  of  nodal  coordinates.  An  example  of  ANCF  finite  elements  is 
the  three-dimensional  gradient  deficient  Euler-Bemoulli  beam  element  shown  in  Fig.  1.  For  this 
element,  the  vector  of  nodal  coordinates  e'A ,  which  contains  position  and  position  gradient 


vectors,  can  be  defined  as  e'A 


driklT 

dxik 


drik2T 

dxik 


(Gerstmayr  and  Shabana,  2006; 


Shabana,  2012),  where  r'k'  and  (dr,kj /dxlk  \  refer  to  the  displacement  and  displacement  gradient 

vectors  at  node  j  of  element  k  .  While  ANCF  finite  elements  were  introduced  for  large 
deformation  problems  of  very  flexible  bodies,  these  elements  can  be  very  effective  in  the 
analysis  of  small  deformation  in  complex  flexible  multibody  systems  containing  a  large  number 
of  joints  and  flexible  links.  As  demonstrated  in  this  paper  using  a  large  scale  tracked  vehicle  that 
has  flexible-link  chains,  ANCF  finite  elements  can  be  used  to  capture  model  details  that  cannot 
be  captured  using  small  deformation  and  rigid  body  formulations. 

A  number  of  approaches  have  been  proposed  in  the  literature  for  developing  tracked 
vehicle  models.  Galaitsis  (1984)  introduced  a  comprehensive  model  incorporating  the  dynamics 
of  independent  rigid  track  shoes.  In  this  model,  the  equations  of  motion  for  the  track  chain  were 
developed  for  a  tracked  vehicle  traveling  on  a  flat  ground.  Bando  et  al  (1991)  developed  a 
procedure  for  the  design  and  analysis  of  small-size  rubber  tracked  bulldozers.  They  presented  a 
computer  simulation  which  was  used  in  the  evaluation  of  the  vehicle  behavior.  Nakanishi  and 
Shabana  (1994)  introduced  a  model  for  the  nonlinear  dynamic  analysis  of  two  dimensional 
tracked  vehicle  systems.  They  used  both  the  augmented  and  the  velocity  transformation 
techniques.  In  their  investigation,  it  was  shown  that  the  singularity  problem,  which  can  be 
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encountered  if  one  set  of  independent  coordinates  is  used  throughout  the  simulation,  can  be 
avoided  by  either  repeatedly  changing  the  set  of  independent  coordinates  or  by  using  a  penalty 
function  approach.  In  this  latter  approach,  a  cut  is  made  at  a  selected  secondary  joint  in  the  track 
in  order  to  form  an  open  loop  structure,  and  the  chain  is  closed  using  a  spring  with  high  stiffness 
coefficient.  Choi  et  al  (1998A;  1998B)  developed  a  three  dimensional  contact  force  model  that 
represents  the  interaction  between  the  track  links  and  the  vehicle  components.  In  their  work, 
track  links  were  connected  using  perfect  pin  joints.  The  highly  nonlinear  joints  constraint 
equations  were  adjoined  to  the  equations  of  motion  using  the  augmented  formulation.  Maqueda 
et  al  (2010)  developed  a  three-dimensional  tracked  vehicle  model  with  rubber  chains  modeled 
using  ANCF  finite  elements  which  allows  for  a  straightforward  implementation  of  general  linear 
and  nonlinear  material  models.  The  cross  section  deformation  results  obtained  using  ANCF  finite 
elements  were  validated  by  comparing  with  the  results  obtained  using  solid  finite  elements  in  the 
case  of  a  simple  tension  test.  In  previous  investigations  (Hamed  et  al,  2011;  Shabana  et  al,  2012), 
an  ANCF  C°/C'  continuity  flexible-link  chain  mesh  was  developed.  This  new  ANCF  mesh  can 
be  used  in  the  efficient  modeling  of  flexible-link  chains.  The  C°/C 1  ANCF  chain  model  leads  to 
linear  connectivity  conditions,  thereby  allowing  for  an  efficient  elimination  of  the  dependent 
joint  variables  at  a  preprocessing  stage.  Using  this  approach,  the  ANCF  mesh  mass  matrix 
remains  constant  in  the  case  of  articulated  systems. 

The  organization  of  this  paper  is  as  follows.  Section  3  discusses  the  problems  of  large 
rotation  and  slope  discontinuities  in  MBS  dynamics.  In  Section  4,  a  brief  review  of  the  recently 
developed  linear  finite  element  ANCF  meshes  that  have  linear  connectivity  conditions  and 
constant  inertia  matrix  is  presented.  This  section  also  proposes  a  new  simple  model  for  the 
sliding  joint.  Section  5  discusses  some  of  the  limitations  of  the  B-spline  representation  when 
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flexible-link  chain  structures  are  considered.  Control  points  used  in  the  B-spline  representation 
may  not  represent  material  points,  leading  to  some  fundamental  issues  that  are  discussed  in 
Section  5.  The  discussion  presented  in  Section  5  sheds  light  on  the  potential  of  using  ANCF  as 
the  basis  for  the  integration  of  computer  aided  design  and  analysis  (I-CAD-A).  Section  6 
provides  a  justification  for  the  use  of  ANCF  finite  elements  to  model  complex  flexible-link 
chains.  Section  7  describes  the  simplified  three-dimensional  Euler-Bernoulli  beam  elastic  force 
model  used  to  demonstrate  the  procedures  and  concepts  introduced  in  this  study.  In  Section  8,  a 
tracked  vehicle  model  is  used  to  obtain  the  numerical  results  that  demonstrate  the  potential  of 
using  the  new  ANCF  meshes  in  capturing  model  details  that  cannot  be  captured  using  existing 
methods.  Summary  and  conclusions  drawn  from  this  study  are  presented  in  Section  9. 

3.  STRUCTURAL  AND  NONSTRUCTURAL  DISCONTINUITIES 

Structural  and  nonstructural  discontinuities  have  been  fundamental  issues  in  the  development  of 
MBS  formulations.  Understanding  these  types  of  discontinuities  will  also  shed  light  on  the 
fundamental  differences  between  the  approaches  currently  being  developed  for  the  integration  of 
computer  aided  design  and  analysis.  Structural  discontinuities  refer  to  the  type  of  joints  that  do 
not  allow  relative  rigid  body  displacements,  while  nonstructural  discontinuities  refer  to  the  types 
of  joints  that  allow  for  relative  rigid  body  displacements.  The  joints  used  in  flexible-link  chains 
are  example  of  joints  that  lead  to  nonstructural  discontinuities. 

Structural  discontinuity  was  the  main  issue  in  the  development  of  the  large  displacement 
small  deformation  FE/FFR  formulation.  In  the  classical  finite  element  formulations,  beams, 
plates,  and  shells,  which  employ  infinitesimal  rotations  as  nodal  coordinates,  are  considered 
nonisoparametric  elements.  Rigid  body  motion  of  such  elements  does  not  lead  to  zero  strains  and 
rigid  body  inertia  cannot,  in  general,  be  exactly  modeled  (Shabana,  2014).  Therefore,  before 
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introducing  the  FE/FFR  formulation  in  the  early  eighties,  an  incremental  approach  was  needed 
when  these  elements  are  used  in  the  analysis  of  large  rotation  problems.  Because  the  incremental 
solution  of  linearized  system  of  equations  eventually  deviates  from  correct  solution,  there  was  a 
need  to  develop  a  new  FE  formulation  that  correctly  describes  rigid  body  motion  when 
nonisoparametric  elements  are  used.  For  this  reason,  the  FFR  formulation  was  introduced.  In  this 
formulation,  a  floating  body  coordinate  system  which  shares  the  large  displacement  and  rotation 
of  the  body  is  introduced.  The  large  displacement  of  the  body  can  be  described  using  the  global 
position  of  the  floating  frame  origin  and  a  set  of  orientation  parameters  that  can  be  used  to  define 
the  FFR  transformation  matrix.  The  use  of  these  independent  reference  coordinate s  is  sufficient 
to  describe  correctly  the  rigid  body  motion.  Nonetheless,  the  introduced  body  coordinate  system, 
if  not  used  with  the  finite  element  method,  allows  only  modeling  simple  geometries  that  do  not 
have  any  slope  discontinuities  as  the  one  shown  in  Fig.  2.  On  the  other  hand,  complex 
geometries  that  consist  of  elements  with  different  orientations  as  shown  in  Fig.  3  can  be  modeled 
using  the  concept  of  the  intermediate  element  coordinate  system.  This  coordinate  system  has  an 
origin  that  is  rigidly  attached  to  the  origin  of  the  body  coordinate  system  and  has  axes  that  are 
initially  parallel  to  the  axes  of  the  finite  element  coordinate  system  as  shown  in  Fig.  4.  Using  this 
coordinate  system,  the  finite  elements  that  have  different  orientations  can  be  connected  rigidly 
using  linear  algebraic  constraint  equations  at  a  preprocessing  stage.  Therefore,  the  FE/FFR 
formulation  provides  a  general  solution  for  accurately  modeling  structural  geometric 
discontinuities  using  linear  connectivity  conditions.  This  leads  to  an  FE  mesh  that  is  locally 
linear  allowing  for  the  use  of  coordinate  reduction  techniques  to  filter  out  high  frequency  modes 
of  vibration. 
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Another  type  of  slope  discontinuity  is  the  nonstructural  discontinuity  which  permits  the 
relative  rigid  body  displacements  between  the  segments  connected.  Example  of  which  is  a  chain 
of  links  connected  by  pin  joints  as  the  one  shown  in  Fig.  5.  In  the  case  of  large  relative  rotation 
between  the  connected  links,  the  dynamics  of  these  chains,  even  for  a  planar  rigid-link  chain,  is 
governed  by  highly  nonlinear  equations  as  the  result  of  the  geometric  nonlinearities  due  to  the 
finite  relative  rotations  (Korkealaakso  et  al,  2009).  For  this  type  of  joints,  the  FFR  or  rigid  body 
formulations  lead  to  a  highly  nonlinear  connectivity  conditions  as  the  result  of  using  orientation 
parameters.  Furthermore,  the  resulting  joint  formulations  do  not  include  any  deformation  modes 
since  the  resulting  joint  degrees  of  freedom  define  only  relative  rigid  body  displacements. 

A  systematic  approach  for  modeling  joints  between  flexible  bodies  using  ANCF  finite 
elements  was  introduced  by  Sugiyama  et  al  (2003).  In  their  work,  tangent  frames  were  used  to 
define  transformation  matrices  at  the  joint  definition  points.  This  approach  leads  to  complex 
nonlinear  constraints  that  are  adjoined  to  the  system  equations  of  motion  using  the  technique  of 
Fagrange  multipliers.  Nonetheless,  as  demonstrated  in  previous  publications,  a  linear  finite 
element  mesh  of  chains  that  consist  of  flexible  links  can  be  developed  in  many  MBS  applications 
(Vallejo  et  al.,  2003;  Hamed  et  al,  2011;  Shabana  et  al,  2012).  This  can  be  achieved  using  ANCF 
finite  elements  which  allow  for  developing  an  FE  mesh  that  includes  one  or  more  chains  with 
flexible  links,  and  each  of  the  links  can  have  independent  relative  finite  rotations.  The 
connectivity  conditions  between  the  ANCF  finite  elements  can  be  described  using  linear 
algebraic  constraint  equations  that  can  be  used  at  a  preprocessing  stage  to  eliminate  redundant 
variables.  The  use  of  this  new  concept  can  lead  to  significant  simplification  of  the  MBS 
algorithm  and  can  eliminate  the  distinction  between  the  augmented  and  recursive  methods.  Both 
the  augmented  and  recursive  methods  will  converge  to  approximately  the  same  CPU  time.  Figure 
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6  shows  two  flow  charts  that  describe  the  steps  of  the  dynamic  analysis  of  two  different  ANCF 
chains  defined  using  the  conventional  nonlinear  constraint  approach  and  the  new  linear  ANCF 
constraint  approach.  It  is  clear  that  using  the  new  ANCF  linear  constraint  approach  is  more 
efficient  since  it  eliminates  many  of  the  time-consuming  steps  of  the  conventional  approach.  For 
instance,  when  the  new  approach  is  used,  there  is  no  need  to  solve  the  constraint  equations  at  the 
position,  velocity  or  acceleration  levels.  Furthermore,  the  dimension  of  the  system 
differential/algebraic  equations  will  be  significantly  reduced  as  a  result  of  the  elimination  of 
dependent  variables  and  the  associated  Lagrange  multipliers. 

4.  ANCF  JOINT  FORMULATIONS  AND  DEFORMATION  MODES 

Using  ANCF  finite  elements,  mechanical  joints  between  finite  elements  can  be  formulated  using 
linear  connectivity  conditions  (Vallejo  et  al.,  2003;  Hamed  et  al,  2011;  Shabana  et  al,  2012).  This 
is  an  important  feature  that  can  be  exploited  to  significantly  reduce  the  number  of  nonlinear 
algebraic  constraint  equations  and  the  associated  Lagrange  multipliers  in  many  MBS 
applications.  Several  flexible  bodies  can  be  modeled  as  one  FE  mesh  consisting  of  elements  that 
have  relative  translations  and  finite  rotations  with  respect  to  each  other.  It  is  demonstrated  in  this 
section,  using  several  examples,  how  the  linear  joint  equations  can  be  developed  and  how  ANCF 
finite  elements  can  also  reduce  the  degree  of  nonlinearity  when  linear  joint  formulations  cannot 
be  obtained  for  some  joints.  The  modes  of  the  deformation  of  the  new  ANCF  joints  are  also 
discussed  in  this  section. 

4.1.  Spherical  Joint 

When  two  bodies  or  finite  elements  are  connected  by  a  spherical  joint,  one  has  C°  continuity 
that  ensures  only  the  continuity  of  the  position  coordinates.  The  two  connected  bodies  can 
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arbitrarily  rotate  with  respect  to  each  other.  For  a  spherical  joint  between  two  finite  elements  i 
and  j  ,  the  algebraic  spherical  joint  constraint  equations  can  be  written  as 

ri,B=rj,A'  (i) 

where  A  and  B  refer  to  the  joint  node  on  the  finite  elements  i  and  j,  respectively.  The  preceding 
constraint  equations  eliminate  the  relative  translation  between  the  two  finite  elements.  However, 
since  ANCF  structural  finite  elements  can  correctly  describe  arbitrary  rigid  body  rotations, 
applying  the  linear  algebraic  equations  of  Eq.  1  at  a  preprocessing  stage  allows  for  developing  a 
FE  mesh  in  which  each  finite  element  can  have  an  independent  rotational  displacement.  Note 
that  the  above  spherical  joint  constraints  are  linear  in  the  ANCF  coordinates.  The  same  constraint 
equations  are  highly  nonlinear  in  the  case  of  the  rigid  body  or  the  FFR  formulation  (Shabana, 
2001). 

4.2.  Pin  Joint 

The  formulation  of  the  three-dimensional  ANCF  linear  pin  joint  will  be  discussed  using  two 
different  ANCF  finite  elements.  The  first  element  is  the  three-dimensional  Euler-Bernoulli  beam 
element.  This  element  has  only  one  gradient  vector  tangent  to  the  element  center  line.  In  order 
for  two  three-dimensional  ANCF  Euler-Bernoulli  finite  elements  to  be  connected  by  a  pin  joint, 
two  sets  of  conditions  must  be  satisfied.  The  first  set  of  conditions  is  the  continuity  of  the 
position  coordinates  at  the  joint  definition  point.  This  set  of  conditions  can  be  imposed  at  a 
preprocessing  stage  by  using  the  same  position  degrees  of  freedom  for  the  two  elements  at  the 
joint  node.  The  second  set  of  conditions  constrain  the  normals  to  the  gradient  vectors  at  the  joint 
point  to  be  parallel  to  each  other.  This  second  set  of  conditions  can  also  be  introduced  at  a 
preprocessing  stage  by  applying  continuity  on  the  component  of  the  gradient  vector  in  the 
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direction  of  the  joint  axis.  The  constraint  equations  of  a  pin  joint  between  two  ANCF  Euler- 
Bemoulli  beam  elements  i  and  j  can  thus  be  written  as  follows: 


where  a  =  1,  2,  or  3  denotes  the  direction  of  the  joint  axis.  The  preceding  equations  represent  a 
set  of  four  scalar  equations.  The  second  equation  in  Eq.  2  is  a  scalar  equation  that  equates  the 
a  th  component  of  element  i  gradient  vector  along  the  joint  axis  to  the  corresponding  a  th 
component  of  element  j  gradient  vector.  For  example,  if  the  joint  axis  is  the  element  y  axis, 


then  in  this  case  a  =  2.  Note  that  the  second  set  of  algebraic  equations  can  be  formulated  at  a 
preprocessing  step  and  can  be  used  to  eliminate  dependent  variables  before  the  start  of  the 
dynamic  simulation.  Using  the  preceding  equations,  one  can  develop  one  linear  FE  mesh  for 
flexible-link  mechanisms  in  which  the  links  can  have  arbitrarily  large  relative  rotations  and  large 
deformations  without  the  need  to  include  nonlinear  algebraic  constraint  equations  in  the 
formulation  since  these  conditions  can  be  eliminated  as  a  preprocessing  stage.  Note  that  for  the 
gradient  deficient  Euler-Bemoulli  beam  element,  the  joint  axis  needs  to  be  one  of  the  finite 
element  axes  in  order  to  obtain  the  linear  relationships  of  Eq.  2.  This  restriction,  which  can  be 
relaxed  when  using  fully  parameterized  ANCF  finite  elements,  does  not  represent  a  serious 
limitation  when  chain  models  are  developed. 

Another  example  of  ANCF  finite  elements  is  the  three-dimensional  shear  deformable 
beam  element.  A  node  on  this  fully  parameterized  ANCF  element  has  12  degrees  of  freedom. 
These  degrees  of  freedom  include  three  rigid  body  translations,  three  rotations,  and  six 
deformation  modes.  One  may  introduce  a  pin  joint  formulation  that  uses  the  following  six  scalar 
equations  defined  at  the  joint  node  between  elements  i  and  j  (Shabana  et  al,  2012). 
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„i.B 

r  =  r 


dr‘’B  _  drjA 
dp  ~  dp 


(3) 


In  the  above  equations,  P  is  the  coordinate  that  defines  the  joint  axis.  It  can  be  either  x,  y,  z,  or 
any  other  coordinate  line  since  the  use  of  ANCF  finite  elements  allows  for  a  straight  forward 
application  of  the  gradient  tensor  transformation.  The  six  scalar  equations  represented  by  Eq.  3 
define  a  pin  joint  with  six  degrees  of  freedom:  one  rotation  and  five  deformation  degrees  of 
freedom.  Consequently,  the  joint  leads  to  C1  continuity  for  the  gradient  vector  along  the 
coordinate  line  p  and  C°  continuity  for  the  remaining  gradient  vectors. 


4.3.  Sliding  Joint 

Sliding  joints  between  flexible  bodies  are  examples  of  variable  boundary  conditions  joints. 
Developing  the  kinematic  equations  of  these  joints  using  the  FFR  formulation  that  employs 
linear  modes  can  be  difficult  because  the  change  of  the  boundary  conditions.  On  the  other  hand, 
these  joints  can  be  systematically  formulated  using  ANCF  finite  elements  as  explained  in  the 
literature  (Sugiyama  et  al,  2003).  Using  ANCF  finite  elements,  a  procedure  was  introduced  to 
model  the  sliding  joints  using  local  frames  defined  at  the  joint  definition  point.  A  tangent  frame 
can  be  used  to  define  a  joint  local  coordinate  system  whose  transformation  matrix  at  the  joint 
definition  point  can  be  evaluated.  In  order  to  define  the  time  derivative  of  the  constraint 
equations,  the  time  derivatives  of  the  columns  of  the  transformation  matrices  are  written  as 
function  of  the  generalized  coordinates  and  velocities.  This  procedure  leads  to  complex 
constraint  equations  that  are  difficult  to  implement.  The  same  joint  constraints  can  be  formulated 
using  a  simple  and  straightforward  approach  by  introducing  two  nongeneralized  coordinates. 
The  first  nongeneralized  coordinate  u  defines  the  local  position  of  the  joint  definition  point  on 
the  base  element  as  shown  in  Fig.  7,  while  the  second  one  m  defines  the  ratio  between  lengths  of 
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the  two  element  position  gradient  vectors  in  the  direction  of  the  sliding  axis.  The  following 
equations  can  be  used  to  define  the  constraint  equations  of  the  sliding  joint: 


C(e'  ,e7 , u' )  =  r'(e',u‘)  -  r 7 (e7 )  =  0,  =  m^— 

p  p  DC  dy 


In  the  above  equations,  £  and  fare  the  coordinate  lines  in  the  direction  of  the  joint  axis  on  each 


element,  and,  e'  and  e7  are  the  nodal  coordinates  of  the  two  finite  elements.  These  constraint 
equations  can  be  applied  to  constrain  4  of  the  joint  degrees  of  freedom:  2  rigid  rotations,  and  2 
rigid  translations.  The  sliding  is  ensured  by  requiring  continuity  of  the  gradient  along  the  sliding 
axis  and  preventing  the  joint  definition  point  from  moving  normal  to  the  direction  of  the  sliding 
axis.  This  guarantees  the  freedom  of  sliding  without  the  need  to  define  local  frames  which  lead 
to  complex  constraint  equations  required  to  ensure  the  parallelism  of  the  sliding  axes. 

4.4.  ANCF  Joint  Deformation  Modes 

In  rigid  body  and  small  deformation  FFR  MBS  algorithms,  joint  formulations  allow  for  only 
relative  rigid  body  degrees  of  freedom  that  describe  the  relative  displacements  between  two  rigid 
triads  attached  to  the  two  bodies  connected  by  the  joint.  The  ANCF  joint  formulations,  on  the 
other  hand,  allow  for  deformation  modes  at  the  joint  definition  point.  This  leads  to  a  new  class  of 
ideal  compliant  joints  defined  using  algebraic  equations.  The  gradient  vectors  of  the  two 
elements  can  change  their  magnitude  and  orientation  allowing  for  stretch  and  shear  deformations 
at  the  joint  node.  Depending  on  the  element  type  and  the  number  of  algebraic  equations  used  in 
the  formulation  of  the  joint,  different  deformation  modes  can  be  captured.  As  an  example,  the 
gradient  deficient  Euler-Bernoulli  beam  element  allows  for  only  the  axial,  bending,  and  twist. 
Therefore,  such  an  element  does  not  capture  shear  and  cross  section  deformations.  On  the  other 
hand,  the  fully  parameterized  beam  element  which  has  a  complete  set  of  gradient  vectors  can  be 
used  to  capture  the  change  in  the  cross  section  due  to  deformation  at  the  joint  node.  For  the  fully 
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parameterized  beam  element,  the  modes  of  deformations  at  the  joint  nodes  can  be  examined 


using  the  Green-Lagrange  strain  tensor  £  which  can  be  written  as  £ 


J -ij,  where  J  is  the 


matrix  of  position  vector  gradients,  and  I  is  the  identity  matrix.  Because  the  strain  tensor  £  is 
symmetric,  one  can  define  the  strain  vector  £r  =  [■£„  syy  szz  cxy  syz  c7X  ] T ,  where  sxx ,  syy  and  szz 
are  the  normal  strains,  and sxy ,  s  and  s7X  are  the  shear  strains.  The  strain  components  are 

continuous  inside  the  element  and  some  of  the  strain  components  can  be  discontinuous  at  the 
joint  nodes  because  of  the  discontinuity  of  position  vector  gradients  at  these  nodes.  Figure  8 
shows  some  of  the  deformation  modes  of  the  three-dimensional  fully  parameterized  beam 
element.  These  modes  are  the  tensile  mode,  the  bending  mode,  and  one  of  the  cross  section 
deformation  modes.  The  tensile  mode  is  a  result  of  changing  the  magnitude  of  the  gradient  rv  at 
the  joint  point.  The  bending  mode  is  due  to  the  change  of  the  orientation  of  the  gradient  vector 
rvwith  respect  to  the  line  connecting  the  two  nodes  of  the  element.  Finally,  the  cross  section 
deformation  mode  defines  the  change  of  the  cross  section  dimensions  which  is  the  result  of  the 
change  in  the  magnitude  of  the  gradient  vectors  along  the  cross  section  fibers.  Furthermore,  the 

second  condition  of  Eq.  3,  (dr'^/sp)  =  (or/,/'/oP) ,  allows  for  the  stretch  of  the  joint  pin. 


5.  COMPARISON  WITH  THE  ISOGEOMETRIC  APPROACH 

Computational  geometry  methods  such  as  B-spline  and  NURBS  (Non-Uniform  Rational  B- 
Splines)  use  control  points  as  coordinates.  The  control  points  do  not  necessarily  represent 
material  points.  This  is  one  of  the  fundamental  differences  between  computational  geometry 
methods  and  the  finite  element  method.  The  isogeometric  analysis  is  a  concept  aims  at  unifying 
the  two  fields  of  computer  aided  design  (CAD)  and  finite  element  analysis  (FEA)  by  using  the 
same  basis  for  the  geometry  description  and  analysis  (Hughes  et  al,  2005).  This  concept,  which 
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has  been  used  in  computer  science  for  many  years  to  build  physics  based  models,  has  recently 
received  attention  from  the  mechanics  community.  The  increasing  interest  in  isogeometric 
analysis  is  motivated  by  the  desire  to  deal  with  the  difficulties  and  cost  associated  with  the 
conversion  of  solid  models  to  FE  analysis  models.  This  conversion  process  leads  to  distortion  of 
the  geometry  of  the  analysis  model  because  existing  structural  elements  such  as  beams,  plates 
and  shells  cannot  describe  correctly  arbitrary  geometry. 

While  the  use  of  the  isogeometric  analysis  can  be  considered  as  a  step  in  the  right 
direction,  the  use  of  computational  geometry  methods  as  analysis  can  lead  to  fundamental  and 
implementation  problems  which  can  be  avoided  by  using  ANCF  geometry.  In  this  section,  some 
of  the  limitations  of  computational  geometry  methods  when  used  as  analysis  tools  are  discussed. 
The  discussion  in  this  section  will  focus  on  the  limitations  associated  with  the  joint  formulations 
which  is  the  subject  of  this  paper. 

5.1  B-Spline  Surfaces 

In  this  subsection,  the  B-spline  representation  is  briefly  reviewed  in  order  to  provide  the 
background  required  for  the  discussion  that  follows  in  subsequent  subsections.  B-spline  surfaces 
can  be  defined  in  the  following  parametric  form  (Piegl  and  Tiller,  1997): 

n  m 

(5) 

i= 0  7=0 

where  u  and  v  are  the  parameters,  Nip{u)  and  N  j  q(v)  are  B-spline  basis  functions  of  degree  p 
and  q,  respectively,  and  P  .  are  a  set  of  (n  +  1 )  x  (rn  +  I )  bidirectional  control  points.  The  B- 
spline  basis  functions  Nj  p(u)  can  be  defined  as 
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Ni0(u)  = 


1  if  Uj  <  u  <  i 
0  otherwise 


H+i 


\t  /  \  U  U:  M i+p+\ 

NiJ  u )  = - Ni,p-l(  u>+ - 


-N 


ui+p-u, : 


ui+p+ 1  u 


i+l,p-l 


(u) 


i+ 1 


(6) 


where  ut  ,  i  =  0,1,2,. .. .,n  + p  +  1  are  called  the  knots,  and  w,.  <ui+l.  The  vector 

U  =  {uQ  w,  .  un .  p .  \ }  is  called  the  knot  vector.  Similar  definitions  can  be  introduced  for 

A if  .  (v)  with  another  knot  vector  V  =  { v0  v1  .  } .  The  knots  of  B-spline  surfaces  do 

not  have  to  be  distinct;  distinct  knots  are  called  breakpoints  and  define  surface  segments  with 
non-zero  dimensions.  The  number  of  the  non-distinct  knots  in  U  and  V  at  a  point  is  referred  to 
as  the  knot  multiplicity  associated  with  the  parameters  u  and  v ,  respectively.  The  degree  of 
continuity  at  a  breakpoint  is  controlled  by  the  multiplicity  of  the  knots  at  this  point.  A  decrease  in 
the  knot  multiplicity  by  one  leads  to  an  increase  in  the  degree  of  continuity  by  one.  This  amounts 
to  automatic  elimination  of  one  control  point.  The  concept  of  the  knot  multiplicity  can  lead  to 
fundamental  implementation  problems  when  ideal  compliant  joints  are  used  in  MBS  algorithms 
as  discussed  later  in  this  paper. 

5.2  Basis  Functions 

ANCF  geometry  does  not  impose  restrictions  on  the  number  of  basis  functions  used  in  the 
interpolating  polynomials.  In  contrast,  B-spline  rigid  recurrence  formula  restricts  the  number  of 
basis  functions  depending  on  the  degree  of  the  interpolating  functions;  there  is  a  relationship 
between  the  polynomial  degree,  the  number  of  knots,  and  the  number  of  control  points.  If  r  + 1  is 
the  number  of  knots  in  U  and  s  + 1  is  the  number  of  knots  in  V  ,  then  in  B-spline  geometry,  one 
must  have 

r  =  n  +  p  +  l,  s  =  m  +  q  +  l  (7) 

Using  Eq.  7,  one  can  show  that  if  a  cubic  interpolation  is  used  for  both  u  and  v ,  as  in  the  case  of 
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a  thin  plate,  the  B-spline  representation  requires  16  control  points  regardless  of  whether  the 
shape  of  deformation  of  the  plate  is  simple  or  complex  (Shabana  et  al,  2012).  This  can  be  a 
disadvantage  in  the  analysis  since  such  a  geometric  representation  can  unnecessarily  increase  the 
dimension  of  the  analysis  model  and  leads  to  a  loss  of  the  flexibility  offered  by  the  FE  method. 
As  the  degree  of  the  polynomial  interpolation  increases,  this  drawback  of  B-spline  becomes  more 
serious. 

5.3  Knot  Multiplicity 

In  computational  geometry  methods,  a  decrease  in  the  knot  multiplicity  by  one  leads  to 
elimination  of  one  control  point  and  to  an  increase  in  the  degree  of  continuity  by  one.  The  knot 
multiplicity  cannot  be  used  to  impose  continuity  condition  on  only  one  element  of  a  position 
vector  or  a  gradient  vector.  In  three-dimensional  analysis,  a  minimum  of  three  scalar  equations 
are  automatically  imposed  when  the  knot  multiplicity  is  decreased  by  one.  Some  ideal  compliant 
joints,  as  discussed  in  this  paper,  may  require  imposing  only  one  or  two  algebraic  equations  to 
obtain  a  certain  degree  of  continuity  for  one  or  two  elements  of  a  position  or  gradient  vector. 
Such  ideal  compliant  joint  constraints  cannot  be  formulated  using  the  rigid  B-spline  and  NURBS 
knot  multiplicity  concept. 

5.4  Structural  Discontinuities 

Structural  discontinuities,  such  as  T-,  L-,  and  V-sections,  which  can  be  thought  of  as  a  sudden 
change  in  geometry,  represent  discontinuities  in  gradients  as  the  result  of  the  discontinuities  of 
the  coordinate  lines  (directions).  Modeling  this  type  of  discontinuity,  which  characterizes  both 
mechanical  and  structural  systems,  requires  the  use  of  gradient  transformation  which  is 
fundamental  in  the  implementation  of  the  ANCF  geometric  representation  (Hamed  et  al,  2011). 
Using  C°  continuity  in  B-spline  analysis  is  equivalent  to  having  a  pin  joint  in  the  planar  analysis 
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and  to  a  spherical  joint  in  the  spatial  analysis.  The  case  of  structural  discontinuity,  which  is  also 
of  the  C°  type,  requires  additional  algebraic  equations  in  order  to  eliminate  the  relative  rotation 
at  the  joint  definition  point.  Therefore,  there  are  important  geometric  configurations  where  the 
degree  of  continuity  does  not  change  by  eliminating  degrees  of  freedom  or  equivalently  control 
points.  These  important  geometric  configurations  cannot  be  covered  by  the  rigid  structure  of  the 
B-spline  and  NURBS  which  relates  the  number  of  degrees  of  freedom  or  the  control  points  to  the 
degree  of  continuity  (knot  multiplicity). 

5.5  Sliding  Joints 

Another  important  geometric  scenario  associated  with  the  ideal  compliant  constraints  is  the 
sliding  joint.  This  type  of  joint  requires  continuity  of  the  gradients,  but  not  continuity  of  all 
position  coordinates.  In  the  B-spline  representation,  C1  requires  first  C°  continuity.  As 
previously  discussed,  ANCF  representation  does  not  require  this  sequence,  and  for  this  reason, 
such  an  ANCF  representation  can  be  used  to  model  the  sliding  joint  in  a  straightforward  manner. 
Modeling  the  same  joint  using  B-spline  geometry  is  not  as  straight  forward:  at  a  given 
breakpoint,  the  multiplicity  associated  with  a  knot  can  be  used  to  change  the  continuity  of  the 
derivatives.  For  example,  in  the  case  of  cubic  interpolation,  C°,C\  and  C2  conditions 
correspond,  respectively,  to  knot  multiplicity  of  three,  two,  and  one.  Applying  higher  degree  of 
continuity  implies  all  lower  degree  continuity.  In  other  words,  a  knot  multiplicity  of  three  which 
corresponds  to  C2  implies  C°and  C1  continuity.  Similarly,  applying  C1  implies  C°  continuity. 
Therefore,  B-spline  structure  does  not  allow  applying  continuity  on  position  gradients  without 
position  continuity.  Therefore,  two  sliding  bodies  must  be  modeled  as  two  separate  meshes  and 
additional  algebraic  constraint  equations  must  be  added  to  guarantee  the  sliding  motion. 
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6.  USE  OF  ANCF  FOR  SMAFF  DEFORMATION  PROBFEMS 

Chains  are  used  in  many  important  MBS  applications,  and  therefore,  it  is  important  to  be  able  to 
develop  efficient  and  accurate  chain  models,  particularly  when  the  link  flexibility  is  considered. 
For  example,  chains  are  critical  in  defining  the  mobility  of  tracked  vehicles  that  consist  of  a  large 
number  of  components  interconnected  by  joints  and  force  elements.  The  dynamics  of  such 
systems  are  governed  by  complex  relationships  as  result  of  the  relative  motions  and  joint  forces. 
Furthermore,  considering  the  flexibility  of  the  track  links  increases  the  system  complexity,  model 
dimension,  and  degree  of  nonlinearity.  For  this  reason,  it  is  necessary  to  carefully  select  the 
procedure  used  to  model  the  flexible  links  and  the  joints  connecting  them  in  a  chain  structure. 

As  previously  discussed,  an  efficient  formulation  that  can  be  used  to  model  complex 
chain  systems  is  the  ANCF  which  allows  for  the  elimination  of  the  joint  constraint  equations  at  a 
preprocessing  stage,  allowing  for  the  elimination  of  the  joint  dependent  variables  before  the 
dynamic  simulation  starts  and  leading  to  a  chain  mesh  that  has  a  constant  inertia  matrix  and  zero 
Coriolis  and  centrifugal  forces.  Furthermore,  the  ANCF  geometry  description,  in  contrast  to 
conventional  finite  elements,  is  consistent  with  the  descriptions  used  in  computational  geometry 
methods  (Lan  and  Shabana,  2010;  Mikkola  et  al,  2013). 
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While  ANCF  allows  for  the  use  of  a  general  continuum  mechanics  approach,  special 
elastic  force  models  can  be  derived  and  implemented  in  ANCF/MBS  algorithms  for  the  solution 
of  specific  applications.  In  the  case  of  the  tracked  vehicle  model  considered  in  subsequent 
sections,  the  deformation  of  the  track  links  is  small.  Nonetheless,  the  development  of  FFR  model 
to  study  the  flexible  link  chain  dynamics  is  not  recommended  because  of  the  following  specific 
reasons: 

1.  The  FFR  formulation  leads  to  a  highly  nonlinear  mass  matrix  and  complex  non-zero 
Coriolis  and  centrifugal  forces.  On  the  other  hand,  ANCF  leads  to  a  constant  inertia 
matrix  and  zero  Coriolis  and  centrifugal  forces.  ANCF  Cholesky  coordinates  can  be  used 
to  obtain  an  identity  generalized  inertia  matrix,  leading  to  an  optimum  sparse  matrix 
structure. 

2.  The  FFR  formulation  leads  to  highly  nonlinear  constraint  equations  for  the  joints  between 
the  chain  links.  These  constraint  equations  must  be  satisfied  at  the  position,  velocity,  and 
acceleration  levels,  and  as  discussed  before,  contribute  to  increasing  the  problem 
dimensionality  and  the  array  storage  needed.  In  the  case  of  flexible-link  chains,  the 
overhead  can  be  significant.  On  the  other  hand,  ANCF  finite  elements  allow  for  the 
formulation  of  linear  connectivity  conditions,  and  therefore,  the  dependent  variables  can 
be  eliminated  systematically  at  a  preprocessing  stage.  The  use  of  the  new  ANCF  meshes 
can  lead  to  significant  simplification  of  the  solution  algorithm  and  significant  reductions 
in  the  array  storage  as  previously  discussed. 

3.  In  both  rigid  and  FFR  chain  models,  the  use  of  closed  chains  can  lead  to  singular 
configurations  if  one  set  of  independent  coordinates  is  used.  This  is  not  the  case  when  the 
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new  ANCF  models  are  used.  No  singularities  are  encountered  because  the  joint  constraint 
equations  and  the  joint  dependent  variables  are  eliminated  at  a  preprocessing  stage. 

4.  In  the  case  of  simple  geometry,  one  chain  link  can  be  modeled  using  one  finite  element. 
Therefore,  the  FFR  coordinate  reduction  methods  will  not  be  advantageous,  particularly 
if  the  high  frequency  axial  modes  need  to  be  included  in  the  model  in  order  to  capture  the 
link  tension  forces.  Furthermore,  in  tracked  vehicle  applications,  the  contact  forces  have 
very  high  frequency  contents  that  can  be  comparable  to  the  frequencies  that  results  from 
the  elastic  forces. 

Because  of  these  reasons,  the  large  deformation  ANCF  finite  elements  can  be  effectively  used  to 
efficiently  solve  small  deformation  problems.  ANCF  finite  elements  allow  for  implementing 
small  deformation  elastic  forces  models  as  described  in  the  following  section. 

7.  SMALL  DEFORMATION  ANCF  ELASTIC  FORCE  MODEL 

In  this  section,  a  simple  force  model,  which  can  be  used  in  the  analysis  of  small  deformation 
based  on  a  three-dimensional  ANCF  Euler-Bemoulli  beam,  is  introduced.  The  developed  force 
model  accounts  for  the  coupling  between  bending  and  axial  deformations.  The  model  can  be 
considered  as  a  generalization  of  the  two-dimensional  force  model  developed  by  Berzeri  et  al 
(2001).  Despite  the  simplicity  of  the  proposed  force  model,  it  accounts  for  the  elastic 
nonlinearity  in  the  strain-displacement  relationship. 

7.1  Generalized  Elastic  Forces 

The  strain  energy  of  the  Euler-Bernoulli  beam  element  can  be  written  as  a  function  of 
longitudinal  strain  ex  and  curvature  k  as 
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U  =  -\[EAex2  +  El ic2]dx  . 


(8) 


In  the  preceding  equation,  x  is  the  longitudinal  coordinate,  l  is  the  length  of  the  element,  E  is 
the  modulus  of  elasticity,  A  is  the  cross  section  area,  and  I  is  the  second  moment  of  area.  The 
expression  for  the  vector  of  elastic  forces  can  be  determined  by  differentiating  the  strain  energy 
with  respect  to  the  ANCF  coordinates  as 
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In  the  computer  implementation,  it  is  possible  to  derive  the  expressions  for  the  derivatives  in  Eq. 
9  by  means  of  a  symbolic  computer  program,  or  to  derive  them  analytically  in  order  to  minimize 
the  number  of  mathematical  operations.  Closed  form  expressions  for  the  derivatives  of  both  the 
axial  strain  and  the  curvature  can  be  written  as  (Gerstmayr  and  Shabana,  2006) 
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The  computation  of  these  complex  expressions  can  be  time  consuming,  especially,  for  models 
containing  a  large  number  of  elements.  Berzeri  et  al  (2001)  introduced  a  simple  model  for  the 
elastic  forces  of  the  planar  Euler-Bemoulli  beam  element.  Their  model  can  be  extended  to  the 
spatial  Euler-Bernoulli  beam  element  assuming  small  longitudinal  deformation,  which  is  the  case 
for  the  steel  track  links  used  in  the  chains  of  the  tracked  vehicle  model  considered  in  this 
investigation.  The  expression  for  the  elastic  forces  due  to  longitudinal  forces  can  be  written  as 
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Q e>/  =  j  EAsx  ^dx  =  J  EAexS'xSxedx . 

A  A 


(11) 


By  assuming  small  variation  of  the  strain  over  the  length  of  the  beam,  the  strain  ex  can  be 

factored  out  from  the  expression  of  the  elastic  forces  , allowing  writing  the  preceding  equation  as 

/ 

Qe  J  =  EAIX\ SrxS  Xdxe  =  KeJe ,  (12) 


where  Kt>/  =  EAsx^ STxSxdx .  The  average  axial  strain  can  be  simply  calculated 


as 


s x=(d-l)/l  where  d  is  the  current  length  of  the  beam. 

The  expression  of  the  elastic  forces  due  to  bending  can  be  written  as 


Qe,r  =\EIk^(1x. 


de 


(13) 


In  the  case  of  small  longitudinal  deformation,  the  curvature  can  be  defined  as  |rvt  | .  Hence,  the 
transverse  elastic  forces  can  be  written  as 


Q«,r  =  E1\ slAAe  =  K  ,,e , 


(14) 
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where  Ke ,  =  El \sTxxSxxdx . 


Finally,  the  elastic  force  vector  can  be  written  as 

Q?  =  (K<-,/  +K<..r)e' 


(15) 


where  Ke/  and  K(,  t  arc  constant  stiffness  matrices  associated  with  longitudinal  and  bending 

elastic  forces,  respectively. 

7.2  Generalized  Damping  Forces 
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The  internal  damping  characteristics  of  most  known  materials  lead  to  energy  dissipation  when 
the  materials  are  subjected  to  initial  disturbances  or  shock  loadings.  Therefore,  developing 
realistic  flexible  MBS  models  requires  the  inclusion  of  the  effect  of  internal  damping.  One  of  the 
advantages  of  ANCF  finite  elements  is  the  ability  of  capturing  modes  of  deformations  that 
cannot  be  captured  using  existing  structural  finite  elements  (Mohamed  and  Shabana,  2011).  In 
some  applications,  some  of  these  ANCF  modes,  which  can  have  high  frequencies,  may  not  have 
significant  effect  on  the  solution.  In  such  cases,  explicit  numerical  integration  schemes  are  forced 
to  select  very  small  time  steps  to  capture  the  high  frequencies.  Therefore,  for  such  applications,  it 
is  desirable  to  damp  out  these  high  frequency  modes  by  using  an  appropriate  visco-elastic 
constitutive  model.  This  section  introduces  a  suitable  damping  model  that  can  be  used  with  the 
simple  elastic  force  model  introduced  in  the  previous  section.  The  damping  model  introduced  is 
also  based  on  the  assumption  of  small  deformation. 

The  proposed  visco-elastic  model  applies  damping  to  the  deformation  associated  with  the 
gradient  vector  rt  and  the  curvature  vector  rvs.  where  s  is  the  arc  length.  In  order  to  properly 
define  the  deformation  and  avoid  damping  out  rigid  body  displacement  modes,  a  tangent  frame 
can  be  used  to  define  an  orthogonal  matrix  using  the  vectors  rx  and  rss  as 
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.  The  power  dissipated  can  be  written  as  a  function  of  the  damping 


stress  a  d  and  the  strain  rate  £  as  (Vallejo  et  al,  2005) 

^=(KWV-  (16) 

1  V 

Following  a  procedure  similar  to  the  one  given  in  (Mohamed  and  Shabana,  2011;  Mohamed, 
2011),  the  power  dissipated  can  be  rewritten  in  the  case  of  the  Euler-Bernoulli  beam  element  as 
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where  ex  and  k  are  the  longitudinal  strain  rate  and  curvature  rate.  By  differentiating  the 

dissipated  energy  with  respect  to  the  absolute  nodal  velocities,  the  generalized  damping  force  can 
be  written  as 


.  de  de 
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where  vl  and  v2  are  two  damping  constants. 

Using  the  assumption  of  small  variation  of  the  strain  over  the  length,  both  the  longitudinal  strain 
and  strain  rate  can  be  assumed  constant  through  the  length.  Hence,  the  longitudinal  damping 
force  can  be  written  as 
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where,  Kdl  =  A(v1eJ(l  +  2ex)) JStTS  Xdx.  Using  the  same  assumption,  the  curvature  and 

o 

curvature  rate  can  be  assumed  constant  through  the  length.  Therefore,  the  transverse  damping 
force  can  be  written  as 
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where.  KdJ  =  I  (v2k /k^SxxTSxxedx .  Then,  the  generalized  internal  damping  force  vector  can  be 
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written  as 

Q,7  =(Krf/(e,e)  +  K^(e,e))e 


(21) 
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These  generalized  damping  forces  can  be  introduced  to  the  ANCF  element  equations.  The 
generalized  damping  forces  associated  with  the  ANCF  Cholesky  coordinates  can  also  be 
obtained  using  the  Cholesky  transformation. 

8.  NUMERICAL  RESULTS 

A  flexible-link  chain  three-dimensional  tracked  vehicle  example  is  used  in  this  section  to 
demonstrate  the  implementation  of  the  concepts  discussed  in  this  paper.  This  example  represents 
a  complex  system  with  large  number  of  bodies  connected  by  large  number  of  revolute  (pin) 
joints.  Furthermore,  the  mathematical  models  of  such  vehicles  represent  a  stiff  numerical 
integration  problem  due  to  high  frequencies  of  the  contact  forces.  Such  vehicle  models  can  be 
used  to  test  the  robustness  of  the  MBS  algorithm  developed. 

8.1  Model  Description 

The  tracked  vehicle  used  is  an  armored  personnel  carrier  that  consists  of  a  chassis  and  two  track 
chains.  Each  track  system  consists  of  an  idler,  one  sprocket,  5  road-wheels,  and  64  track  links.  A 
schematic  drawing  of  the  selected  tracked  vehicle  is  shown  in  Fig.  9.  Table  1  shows  the  inertia 
properties  for  all  the  tracked  vehicle  model  components  used  in  this  simulation.  The  vehicle  has 
a  suspension  system  that  consists  of  a  road  arm,  placed  between  each  road  wheel  and  the  chassis, 
and  a  shock  absorber  that  consists  of  a  spring  and  a  damper  connecting  each  road  arm  to  the 
chassis.  In  order  to  control  the  tension  of  the  two  tracks,  each  idler  is  connected  to  the  chassis 
through  a  tensioner.  Table  2  shows  the  stiffness  coefficients  k  and  the  damping  coefficients  c  of 
the  contact  models,  the  suspension  system,  and  the  tensioner;  as  well  as  the  friction 
coefficient  //  of  the  contact  models.  Ideal  pin  joints  are  used  to  connect  the  road  arms  and  the 
sprockets  to  the  chassis,  and  the  road-wheels  to  the  road  arms.  The  normal  contact  forces  are 
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calculated  based  on  a  linear  dissipative  contact  force  model  (Choi  et  al,  1998;  Maqueda  et  al, 
2010;  Dufva  et  al,  2007). 

The  models  discussed  in  this  section  are  assumed  to  have  flexible  track  links.  The 
flexible  links  are  defined  as  three-dimensional  steel  Euler-Bemoulli  beam  elements  with  a 
modulus  of  rigidity  of  76.9  GPa,  a  Young’s  modulus  of  200  GPa,  and  a  mass  density  of  7800 
kg/m  .  The  flexible  chain  links  are  connected  to  each  other  using  ANCF  ideal  compliant  pin 
joints,  which  define  a  linear  ANCF  mesh.  As  previously  mentioned,  one  of  the  advantages  of 
using  this  ANCF  definition  of  joints  is  having  linear  joint  constraints  which  allow  for  the 
elimination  of  the  dependent  joint  variables  at  a  preprocessing  stage.  As  a  result,  this  procedure 
avoids  the  computational  cost  required  to  formulate  the  nonlinear  pin  joint  constraint  equations 
at  the  position,  velocity,  and  acceleration  levels.  In  the  case  of  the  tracked  vehicle  used  in  this 
study,  the  proposed  approach  allows  for  the  elimination  of  126  pin  joints  and  630  highly 
nonlinear  constraint  equations  at  a  preprocessing  stage.  In  conventional  MBS  algorithms,  this 
large  number  of  highly  nonlinear  constraint  equations  must  be  satisfied  at  position,  velocity,  and 
acceleration  levels.  If  ANCF  joints  are  used,  all  these  constraints  become  linear  and  can  be 
eliminated  at  a  preprocessing  stage. 

Three  different  tracked  vehicle  models  are  used  to  demonstrate  the  effectiveness  of  using 
the  ANCF  linear  ideal  compliant  joint  model.  In  the  first  vehicle,  the  chain  links  are  modeled  as 
rigid  bodies.  In  order  to  solve  the  singularity  problem  encountered  in  the  case  of  closed  chain 
(Nakanishi  and  Shabana,  1994),  the  first  and  last  links  of  the  chain  are  connected  through  a  very 
stiff  spring.  The  links  of  the  second  and  third  tracked  vehicle  models  are  modeled  as  flexible 
links  connected  using  ANCF  continuum-based  joints.  The  chain  of  the  second  model  has  a 
secondary  joint  by  making  a  cut  and  using  a  very  stiff  spring  to  close  the  chain  as  in  the  case  of 
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the  rigid-link  chains.  In  the  case  of  the  third  model,  on  the  other  hand,  an  ANCF  closed  loop 
chain  is  used.  The  third  model  is  used  to  demonstrate  the  fact  that  using  the  new  linear  ANCF 
joint  model  does  not  lead  to  any  singularity  problems  as  in  the  case  of  conventional  MBS 
nonlinear  joint  algorithms.  An  explicit  Adams -Bashforth  predictor-corrector  method  was  used 
for  the  numerical  solution  of  the  system  differential  equations. 

8.2  Rigid  Body  Motion 

The  angular  velocities  of  the  two  sprockets  of  all  the  three  tracked  vehicle  models  were  specified 
to  follow  the  pattern  shown  in  Fig.  10.  The  response  of  the  chassis  to  this  input  sprocket  velocity 
is  shown  in  Figs.  11-13.  Figure  11  shows  the  forward  velocity  of  the  chassis.  As  expected,  the 
velocity  of  the  chassis  follows  the  same  pattern  as  the  sprocket  velocity.  The  results  of  Fig.  12 
and  13  show  the  forward  and  vertical  positions  of  the  chassis  center  of  mass.  It  can  be  noticed 
from  the  results  presented  in  Fig.  14  that  the  chassis  has  0.03m  downward  displacement  as  a 
consequence  of  an  initial  settlement  of  the  suspension  springs.  These  springs,  in  addition  to  the 
contact  forces,  deform  under  the  effect  of  the  weight  of  the  vehicle  components  until  equilibrium 
positions  are  reached.  From  the  results  presented  in  Figs.  11  to  13,  it  is  clear  that  the  rigid  body 
motion  results  of  the  two  flexible  chain  models  are  in  a  good  agreement  with  the  results  of  the 
rigid  body  model.  Figure  14  shows  the  forward  position  of  an  arbitrarily  selected  point  (end  point 
of  link  31).  At  the  beginning  of  the  simulation,  this  node  is  in  contact  with  the  ground.  After  that, 
the  vehicle  starts  moving,  and  this  point  moves  and  it  comes  into  contact  with  the  idler.  It 
continues  moving  until  it  comes  again  into  contact  with  the  ground.  This  motion  is  repeated  as 
the  vehicle  continues  to  move  forward.  Figure  15  shows  the  motion  trajectory  of  the  same  node 
defined  in  the  chassis  coordinate  system. 

8.3  Forces  and  Deformations 
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In  tracked  vehicle  dynamics,  it  is  important  to  study  the  reaction  forces  at  the  pin  joints.  These 
forces  can  be  used  to  further  analyze  the  strain  and  stress  state  of  the  chain  and  determine 
whether  or  not  the  chain  can  sustain  the  forces  without  failure.  Figures  16  and  17  show  the  joint 
constraint  forces  of  the  first  joint  (marked  in  Fig.  21)  in  the  longitudinal  and  the  vertical 
directions  for  the  rigid  and  ANCF  joint  models.  One  advantage  of  including  the  flexibility  of  the 
track  links  is  to  be  able  to  calculate  online  the  strains  and  stresses.  The  online  computation  of  the 
track  link  deformation  leads  to  a  more  realistic  prediction  of  the  system  dynamics,  this  is  due  to  a 
nature  of  the  coupling  between  different  displacement  modes  and  simulation  parameters.  Figures 
18  and  19  show  the  axial  and  transverse  deformations  of  the  first  link  of  the  right  chain  at  t  =  5s . 
In  these  figures,  the  deformation  is  magnified  by  a  factor  of  10s.  The  deformation  of  this  link 
was  calculated  using  a  local  tangent  frame  that  has  its  longitudinal  axis  parallel  to  the  first  node 
position  gradient  vector.  Figure  20  shows  the  axial  stress  of  the  right  track  chain  at  t-5s,  while 
Fig.  21  shows  the  configuration  of  the  chain  at  t  =  5s .  It  is  clear  from  the  results  of  Figs.  20  and 
21  that  the  chain  span  can  be  divided  into  three  sections.  First  section,  which  is  located  before 
the  sprocket  and  is  directly  affected  by  the  sprocket  forces,  is  the  tight  section  which  has  the 
highest  tension.  The  second  section  is  the  slack  section  located  after  the  sprocket  and  has  the 
lowest  degree  of  tension.  The  third  section  is  a  transition  section  between  the  tight  and  slack 
sections  of  the  chain. 

While  all  the  results  of  rigid  body  motion  of  the  different  models  presented  show  good 
agreement,  the  computational  time  of  the  flexible  model  is  almost  four  times  the  computational 
time  of  the  rigid  tracked  vehicle  model.  This  can  be  justified  because  of  the  high  frequencies 
associated  with  the  track  link  deformation.  The  presence  of  high  frequencies  can  be  confirmed 
by  considering  the  number  of  function  evaluations  of  the  flexible  and  the  rigid  tracked  vehicle 
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models.  The  ratio  between  the  numbers  of  function  evaluations  for  the  rigid  and  flexible  tracked 
vehicle  is  approximately  1 :  8.5 .  Considering  both  the  simulation  time  and  the  number  of  function 
evaluations,  it  can  be  concluded  that  using  the  new  ANCF  linear  mesh  in  modeling  the  flexible 
link  tracked  vehicles  leads  to  significant  computational  saving.  This  computational  saving  is  the 
result  of  the  elimination  of  the  joint  constraints  and  the  corresponding  dependent  coordinates  and 
Lagrange  multipliers.  Furthermore,  numerical  experimentation  showed  that  when  a  low  modulus 
of  elasticity  is  used,  flexible-link  tracked  vehicle  models  become  more  efficient  than  the  rigid- 
link  tracked  vehicle  model  as  the  result  of  using  the  new  ANCF  finite  element  mesh. 

9.  SUMMARY  AND  CONCLUSIONS 

In  this  paper,  a  new  class  of  ideal  compliant  joints  that  account  for  the  distributed  inertia  and 
elasticity  is  introduced  and  used  to  discuss  fundamental  issues  related  to  the  integration  of 
computer  aided  design  and  analysis  (I-CAD-A).  Additional  degrees  of  freedom  are  introduced  in 
order  to  capture  modes  of  deformation  that  cannot  be  captured  using  existing  formulations  that 
employ  rigid  triads  or  unit  vectors  to  formulate  the  joint  equations.  The  ideal  compliant  joints 
developed,  while  more  general  and  capture  more  modes  of  displacements,  can  be,  for  the  most 
part,  formulated  using  linear  algebraic  equations.  These  linear  connectivity  conditions  allow  for 
the  elimination  of  the  dependent  variables  at  a  preprocessing  stage,  thereby  significantly 
reducing  the  problem  dimension  and  array  storage  needed.  Consequently,  the  constraint 
equations  are  automatically  satisfied  at  the  position,  velocity,  and  acceleration  levels.  Therefore, 
when  using  the  proposed  approach  to  model  large  scale  chain  systems,  differences  in 
computational  efficiency  between  the  augmented  formulation  and  the  recursive  methods  are 
eliminated,  and  the  CPU  times  resulting  from  the  use  of  the  two  formulations  become  similar 
regardless  of  the  complexity  of  the  system.  Furthermore,  the  elimination  of  the  joint  constraint 
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equations  and  the  associated  dependent  variables  contribute  to  the  solution  of  a  fundamental 
singularity  problem  encountered  in  the  analysis  of  closed  loop  chains  and  mechanisms,  and  there 
is  no  need  to  repeatedly  change  the  chain  or  mechanism  independent  coordinates. 

Using  the  concept  of  the  ideal  compliant  joints,  the  fundamental  problems  that  can  be 
encountered  when  using  computational  geometry  methods,  can  be  explained.  The  concept  of  the 
knot  multiplicity  used  in  computational  geometry  methods  is  not  suited  for  the  formulation  of 
ideal  compliant  joints  that  may  require  different  degrees  of  continuity  for  different  variables  or 
elements  of  the  gradient  vector  or  require  imposing  constraints  on  the  gradient  vectors  along 
directions  that  differ  from  the  B-spline  and  NURBS  coordinate  lines.  As  explained  in  this  paper, 
this  issue  is  closely  related  to  the  inability  of  B-spline  and  NURBS  to  model  structural 
discontinuities . 

Another  contribution  of  this  paper  is  demonstrating  that  large  deformation  ANCF  finite 
elements  can  be  effective  in  solving  small  deformation  problems.  The  necessity  of  using  the  new 
concepts  developed  in  this  paper  in  modeling  complex  systems  is  demonstrated  using  a  heavily 
constrained  tracked  vehicle  with  flexible  link  chains.  Without  using  the  proposed  approach 
modeling  such  a  complex  system  with  flexible  links  can  be  a  challenging  task.  The  analysis 
presented  in  this  paper  also  demonstrates  that  adding  significant  model  details  does  not 
necessarily  imply  increasing  the  complexity  of  the  MBS  algorithm.  In  fact,  including  significant 
details  such  as  body  flexibility  and  joint  deformations  can  be  accomplished  using  MBS 
algorithms  that  are  much  simpler  than  the  algorithms  currently  used  for  rigid  body  system 
analysis. 
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Table  1.  Tracked  vehicle  components  inertia  properties 


Part 

Mass  (kg) 

Ixx  (kg  m2) 

Iyy  (kg  m2) 

Izz  (kg  m2) 

Ixy,  Iyz,  Ixz  (kg  nr) 

Chassis 

5489.24 

1786.92 

10449.67 

10721.22 

0 

Sprocket 

436.67 

13.87 

12.22 

12.22 

0 

Idler 

429.57 

14.70 

12.55 

12.55 

0 

Road  Wheel 

561.07 

26.06 

19.82 

19.82 

0 

Road  Arm 

75.26 

0.77 

0.37 

0.77 

0 

Track  link 

18.02 

0.04 

0.22 

0.25 

0 

Table  2.  Tracked  vehicle  contact  and  suspension  parameters 


Parameters 

Sprocket-Track 

Contact 

Roller-Track 

Contact 

Ground-Track 

Contact 

Suspension 

Tensioner 

k  (N/m) 

2xl06 

2xl06 

2xl06 

lxlO6 

lxlO6 

c  (N's/m) 

5xl03 

5xl03 

5xl03 

1.5xl04 

1.4xl04 

A 

0.15 

0.10 

0.30 

- 

- 
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Figure  1.  Euler-Bemoulli  beam  element 

Body  coordinate 


Figure  2.  Structure  without  slope  discontinuity 
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Figure.  4  Intermediate  coordinate  system 
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Figure  5.  Eight-link  chain 
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Read  initial  data  including  initial 
estimate  for  the  system  state 

i 

Evaluate  the  AN( 

mai 

- 1 

3F  constant  mass 
trix 

1 - 

(a) 


K 


(b) 


Figure  6.  Algorithm  flow  chart 
(a)Nonlinear  constraints  (b)  New  linear  constraints 
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(a)  Tensile  mode  (b)  Bending  mode  (c)  Cross  section  deformation 
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Figure  9.  Tracked  vehicle  model 


Time  (s) 


Figure  10.  Sprocket  rotational  velocity 
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Figure  11.  Chassis  forward  velocity 

( -  Rigid  model, -  Closed  loop  flexible  model , - Open  loop  flexible  model  ) 


Figure  12.  Forward  displacement  of  the  chassis 
( -  Rigid  model, -  Closed  loop  flexible  model , - Open  loop  flexible  model  ) 
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Figure  13.  Vertical  displacement  of  the  chassis 
( -  Rigid  model, - Closed  loop  flexible  model , -  Open  loop  flexible  model  ) 


Figure  14.  Node-32  forward  position 
( - Rigid  model,  - Flexible  models) 
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Figure  17.  First  joint  vertical  forces 
( - Rigid  model,  -  Flexible  models) 


Y  (m)  0  -0.02  X(m) 


Figure  18.  Axial  deformation  of  first  link  of  the  right  chain  at  t  =  5s  magnified  by  a 

factor  of  105 

( - Initial  configuration, - Deformed  configuration) 
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Figure  19.  Transverse  deformation  of  first  link  of  the  right  chain  at  t  =  5s  magnified  by  a  factor 

of  10' 

( -  Initial  configuration,  - Deformed  configuration) 


x  10'5 


Figure  20.  Right  chain  axial  stress  at  t  =  5s 
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Figure  21.  Right  chain  axial  force  distribution  at  t  =  5s 


